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ABSTRACT 


This thesis investigates the feasibility of micro- 
computer based simulation of scalar wave propagation in var- 
ious media. Models for lossless media and media with a loss 
coefficient which is linear in frequency have been coded om 
FORTRAN and simulated successfully on a commercially avail- 
able micro-computer, with simulation times less than Thess 
minutes. The spatial impulse responses for classical 
problems using square and circular-piston excitation are 
presented graphically, along with new, innovative, spatial 


excitation source shapes. 
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I. INTRODUCTION 


Due to diffraction effects, acoustic imaging and tissue 
characterization techniques require information about the 
insonifying wave at the object's location. The propagation 
of a pulsed ultrasound wave with arbitrary temporal and 
spatial shape is not well understood, and requires 
development of computer-aided analysis and simulation 
techniques. 

A computationally efficient method to model ultrasound 
propagation from a planar source within, a medium has been 
developed by Guyomar and Powers [Refs. 1,2]. The method ‘is 
Bepiicabile to lossless media, media with a linear frequency 
loss coefficient such as biological tissue, and media with a 
quadratic frequency loss coefficient such as liquids and 
gases. Relying on Fast Fourier Transforms (FFTs) instead of 
integral solutions that require geometric interpretation 
makes the method very suitable for computer implementation. 

Models for all three cases have been previously 
developed, coded in FORTRAN, and simulated on an IBM 3033 
mainframe computer [Ref. 3]. The purpose of this thesis was 
to investigate the feasibility of generating micro-computer 
solutions of the models within a realistic amount of time. 
As the average execution time was thirty seconds on the 


mainframe for a 64 x 64 array size, a time limit of thirty 


minutes was set as a realistic goal. While the difference 
between thirty seconds and thirty minutes may appear 
unacceptable, rapid advances in micro-computer technology is 
narrowing this margin appreciably. 

Since the previously written FORTRAN programs were un- 
available, the project required starting from the beginning. 
A first attempt was made using a commercially available pro- 
gram called MATLAB. MATLAB iS an array-oriented program 
containing built-in functions such as Bessel and Fast 
Router Transforms, as well as the ability to program new 
functions. Limitations later discovered in MATLAB required 
starting over, programming the entire project in Microsoft 
FOre ran = Versi One sol. The. result is a stand-alone, 
user-interactive program allowing customized computer runs 
tailored to a specific set of input variables. Only the 
lossless and linear loss coefficient cases have been imple- 
mented, the third case with quadratic loss coefficient is 
left for future expansion. The goal of thirty minutes was 
obtained by code optimization and exploiting array symmetry 
where it existed. 

The program was tested and results verified using clas- 
Sical problems such as a square or circular-shaped piston 
spatial excitation. Once verified, the program was used to 
compute the spatial impulse response of new excitation 
shapes. Examples of these include pyramid-shaped pulsed 


sources and pulsed linear array elements. 
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The problem solved by Guyomar and Powers [Refs. 1,2] 


uses the geometry of Figure l. 


OBSERVATION 
era le 





A _ SOURCE 


Figure 1. Source and Observation Planes. 
Given a separable source of excitation 


v(z,y,0,t) = s(z,y)T(t), (1) 


maeeproblem is to find the acoustic potential ¢(x,y,z,t) at 
an observation point located a distance z above the source 


plane. A rigidly baffled source plane and linear, 


homogeneous media are assumed. It has been shown (Ref. 2] 


that the potential is given by 
$(z,y,z,t) = s(z,y)T(t)2 ye 9(z,y,2,¢). (2) 


The * indicates convolution over the variable appearing 
directly below it and g(x,y,z,t) is the impulse response or 
Green's function that solves the wave equation and meets the 
applicable boundary conditions, as shown in the block 


diagram of Figure 2. 






Propagation 
f 





boundary conditions 


Figure 2. Block diagram of impulse response. 


In the Spatial domain, ¢(x,y,Zz,t) is also given [Refs. 


4-8] as 


o(z,y,2,t) = T(t): A(z, y,2,) (3) 


where h(x,y,z,t) is the "spatial impulse response" and is 


equal to 
h(zyy3z,t) = s(a3u) 5, giz ye) (4) 


where g(x,y,z,t) is the Green's function (or impulse re- 


sponse). In a linear system, the relationship between the 
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spatial impulse response, h(x,y,zZ,t), and the Green's 
function, g(x,y,z,t), is illustrated by the block diagram of 


Figure 3. Using the two-dimensional spatial Fourier trans- 









s(x, y}6 (t) Propagation 
ap 


boundary conditions 


hixy.z,t) = 
slxyhivtalxy.z.t) 


Figure 3. Block diagram:of spatial impulse response. 


form of Equation 4, ¢(x,y,z,t) can be written as 


a 


(x,y, z,t) = T(t)? F-* {89} | (5) 


where the tilde notation indicates the transform of the spa- 
fEeal function. 
In the spatial frequency domain, the transform of the 


Spatial impulse response is 
h(z,y,z,t) =F {39}. (6) 


maem COnvOlutions in Equations 2 and 4 are difficult to 
implement on a computer. The technique presented by 
Equations 3, 5, and 6, and shown graphically in the block 


diagram of Figure 4, reduce two of the convolutions to 


fmetiplication. 


fl 


enw 


From Equation 6 we see that if we interpret S as the 
angular spectrum of s, then @ can be thought to be the 


"propagation transfer function" associated with propagation 






s(x,y) T(t) 





Propagation 
+ 


boundary conditions 


} (x,y,z, t) 

= Tit) ¥ hixy,z.tl 

= slxy)TIt) *** glxy.z,t) 
xyt 


Figure 4. Block diagram of general solution. 


in the medium. This propagation transfer function behaves as 
a time-varying spatial ig ID EKene that increasingly attenuates 
the higher spatial frequencies of the source as time 
increases. 

Three propagation models have been developed, one for 
each type of media. The general technique followed for each 
model begins with the representative wave equation for the 
specific medium and finding the Green's function (or the two 


dimensional spatial transform of the Green's function) which 


solves the wave equation. Using this Green's function or 
its spatial transform, the spatial impulse - response, 
h(x,y,z,t), is found from Equation 6. From this, the acouem 


tic potential ¢4(x,y,z,t) can be calculated for an arbitmamee 


plane in the positive-~-z half space using Equation 5. 


tee 


A. CASE 1: LOSSLESS MEDIA 

The results of the transfer function approach for loss- 
less media follow, as published by Guyomar and Powers [Ref. 
Ve 


The wave equation is the Helmholtz equation 


1 07¢ 
2 ; = 
Oe a (7) 


The Green's function is 


é(ct —R 
g(z,y,2,t) = te (8) 


where 


R= \/r? + y? + 22. (9) 


For this problem the spatial impulse response is 


2s(2, y) 2 5 8[t — (R/e)) 


h = 10 
(2; Y, Z, ) °nrR ( ) 
Taking the spatial transform yields 
oo EA ih (ov c2t? — z?] H(ct — z), (11) 
vie 


where the tilde indicates the two-dimensional spatial 


transform and 


p= fit f. (12) 
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The multiplicative constants can be dropped since the re- 
sults will be normalized to maximum values. The transfer 


function for propagation in lossless media is then 
Gp1 = Jo (pv c*t? — 2?) H(ct — z), (13) 


whemenrd is the Bessel function of the first kind for ona 


zero. For programming considerations the term 
2 = Cte | (14) 


1s calculated as a unit and identified as ZPRIME in both 
this paper and the source code. 

For a given value of z, one can calculate the value of 
Sp1 for each value of time, then take the inverse spatial 
transform of the product. The result is the spatial impulse 


response of Equation 8. 


Be CASE... LOSS COEFFICIENT LINEAR SiNy-REOuUr Se. 

For media which exhibits a loss coefficient that in- 
creases linearly with frequency, the results of Guyomar and 
Powers [{Refs. 1,2] follow. 

The wave equation for media with a linear loss coeffi- 


cient is the telegrapher's equation 


Z 
1 ny ee (15) 


Veo 
? c- Ot? Ot 
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For this case it is the spatial transform of the 
Green's function that is required. It has been found by an 


alternate method [Ref. 9] to be 


Giz f »*) t) = aa 
ep | Vo — (A2c?/4) Vc2t? — z ‘\e e~Ac?t/2 H(ct —z) for p < Ac/2 


(16) 
Jo | Vo — (A2c?/4) Vc?t? — mel e~Ae?t/2 (ct —z) for p > Ac/2, 


where I, is the Modified Bessel function of the fingsct kane 
of order zero. By definition, the propagation transfer 
function is equal to the transform of the Green's function, 


SO 


Gp2(fz, fy, Zz; £) = 
Io Vp? — (A*%c?/4) Vc*t? — z?] e~Ac*t/2 (ct — z) for p < Ac/2 


Jo Rx = (A2e- (4) c-t* — z3] eAe*t/2 (ct — z) for p > Ac/2. 


fit), 


The transform of the spatial impulse response is 
iy fy, 25 t) aoe 
Bleie, allo Vp? — (A2c?/4) Vc? t? — 4] e~Ac*t/2 H(ct —z) for p< Ac/2 


18 
Sl hes taJ0 Vo — (A*c?/4) Vc?t? — Z| e~Ac*t/2 Het — z) for p > Ac/2. = 


The algorithm for finding the spatial impulse response 
1s similar to the lossless case. The transform is found of 
the driving function s(x,y) and multiplied by Tp2 evaluated 
at each spatial frequency. Taking the inverse transform of 


this product yields the required spatial impulse response. 


S 


III. PROGRAM DESCRIPTION 


A. PROGRAMMED IN FORTRAN 

Fortran was the language chosen for the program for 
three reasons: 1) the availability of International 
Mathematical and Statistical Library (IMSL) (Ref. 10] rere 
tines, specifically those used to calculate Bessel functions 
and one- or two-dimensional Fast Fourier Transforms, 2) the 
portability of the source code to other brands of micro- 
computers and compilers, and 3) familiarity of the source 
code for other students and professionals. Microsoft 
Fortran Version 3.31 was the specific brand of compiler 


used. 


B. ADVANTAGE OF SYMMETRY 

Most of the arrays generated by the program are sym- 
metrical. This was used as an advantage in reducing the 
number of calculations required and therefore the execution 
time of the program. As operation of the program is de- 
scribed, the areas where symmetry has been used are identi- 
fied. Figure 5 is provided to aid in these explanations by 
providing a "map" of an array base. Figure 5 shows a 64 x 
64 array base situated on the X,Y plane and divided into 


four quadrants. Note that the quadrants are unequal in 


16 


size, i.e., quadrant II includes rows 1 through 33 and 


columns 1 through 33 for a 33 x 33 "sub-array" while quad- 
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Figure 5. Base array configuration. 


rant I includes rows 1 through 33 and columns 33 through 64 
Poeea 32 X 32 "sub-array". Similarly, when the left side 
(quadrants II and III) is to be calculated as a whole, the 
calculations cover all 64 rows and columns 1 through 33. 
The different sizes are required to ensure that only column 
Bemeoncains “peak" information and also that element (33,33) 


is the physical center of the array. The importance of this 
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and the reason why it is required will be explained later. 
In-line code provides the required "flipping" actions to 
expand one or two quadrants into a full 64 x 64 array of in- 


formation. 


C. BLOCK DIAGRAM 
Figure 6 is a block diagram illustrating the infor- 
mation flow through the program. The program is explained 
block by block in the following text. 
1. User input 
The program starts by obtaining information specific 
for the problem to be run. Using on-screen prompts, the 
user selects a problem name, filter type, excitation shape 
and size, and values Fox Z, RHO, and MAXTIME. If filter Sp2 
1s selected for lossy media, the user is further prompted 
for a value of the loss coefficient A. 
| The problem name assigned must be eight characters 
or less and conform to MS-DOS filename specifications. Once 
entered, the program concatenates three different extensions 
to the filename. This creates the names of the three output 
files generated by the program. 
The first file created is <filename>.IMP, where 
<filename> is the filename entered by the user. The program 
stores the 64 x 64 spatial excitation array, s(x,y), in the 


default directory under this name prior to taking its two 
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dimensional FFT. It is provided should the user want to 


graphically show the spatial excitation shape. 


The second file created is <filename>.TXT. This is 
a printable text file which summarizes the user's inputs for 


a particular problen. 





/FILTER Spt 
“LossLess OSSY (LINEAR) 
IMEDIA : 


L_ 
Mie WA 








Figure 6. Program block diagram. 


The third and last file created is <filename>.DAT. 
[meee tile contains the final 64 x 64 array which is the spa- 


tial impulse response for the given problem. 


All three files are in ASCII text format. They can 
be printed using the DOS PRINT command, put into a word 
processing program for printing, or put into a graphics 
program for conversion into a picture. All figures 
generated for this thesis were produced by translating the 
<filename>.IMP or <filename>.DAT file to MATLAB format using 
the translation program supplied with MATLAB. The converted 
files were then read into MATLAB and plotted using the MESH 
command. As new, improved, graphic routines for micro- 
eee tage are introduced, a future ihprovenone would be 
incorporation of graphic routines into the program. 

As the program is written, fourteen different source 
excitation configurations are allowed, nine of which have 
been implemented. The choice of fourteen different configu- 
rations was arbitrary, providing a good selection of choices 
and ease of future expansion. Each of the remaining five 
configurations can be implemented by adding its name to the 
screen formatting code and writing the FORTRAN code to gen- 
erate the excitation shape. The desired source excitation 
shape is selected by entering the selection number found 
next to the shape's description. 

The variable RHO represents the radial distance in 
the spatial frequency domain as shown in Equation 12. It is 
identified by the Greek character p in Equations 11 through 
13 and 16 through 18 but will be referred to as RHO to match 


the source (code, RHO is used to create a1 x 64 vector 
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called RHO1 which in turn is used to generate a 64 x 64 ar- 
ray called RHO2. The vector RHO1 is unique since it must 
start approximately at -RHO and increment until it is equal 
to 0.0 in RHO1(33). From here it continues to increase un- 
til it equals +RHO in RHO1(64). Since the vector is of even 
PActh, it cannot be formed simply by dividing the entire 
range of RHO (-RHO to +RHO) by 64 and using this value as 
the increment. To do so would create equal values in 
RHO1(32) and RHO1(33) which defeats the requirement that 
RHO1(33) contain singular information, in this case the 
value 0.0. 

The largest value of time, in seconds, for the time 
dimension is peeneceneea by the variable MAXTIME. For each 
problem, time starts ae the instant that the wave front 
first arrives at the observation plane. This is when time = 
Z/C, where C is the sound velocity (1500 meters/second), and 
increments linearly until it equals MAXTIME. To prevent a 
run-time error, the user is prompted for a value of MAXTIME 
which is Brake: than Z/C seconds. 

The height of the observation plane above the X,Y 
meame tS represented by Z (in meters). It is typically a 
small number, several hundredths of aemecen, On can "be set 
moe zero. A value of 0.1 meters (10 cm) was used for 
debugging the program and producing the results found in the 


Numerical Simulations section. 
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2. Generate ZPRIME vector 

At the completion of user input, the program starts 
to generate the required vectors and arrays. The first of 
these, the TIME vector, is initialized and then used to 
generate the ZPRIME vector. Both vectors are 1 x 64 in size 
and represent a series of time steps starting at time = Z/C 
and increasing to MAXTIME. 

The program contains a variable called STEP. Its 
default value of 4 (which can be changed) sets rows 1 


through 4 of the final RESULT array to zero after all 


calculations have finished. This simulates the step 
function H(ct-Zje Since any information contained in the 
first four rows is lost, these values need not be 


calculated. For this reason, TIME starts at TIME(STEP+1), 
or TIME(5), with a value of 2Z2/C seconds, and increments 
linearly.until it reaches MAXTIME in TIME(64). The smooth 
linear progression of time from Z/C to MAXTIME is obtained 
using the following code segment: 

INC = (MAXTIME - TS) /REAL(NROWS-STEP-1) 

DO XX I = (STEP+1) ,NROWS 

TIME(I) = TS + (I-STEP-1) * INC 
CONTINUE 


where TS = Time Start = Z/C, NROWS 


64, and STEP = 4. 

AS an example let MAXTIME = 1.0E-4 seconds and Z = 
0.1 meters as determined by the user. The quantity TS 
equals 6.666667E-5 seconds and this value is stored in 


TMG Sees Each element of TIME 1S increased by INC = 
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5.649718E-7 until it equals MAXTIME in TIME(64). The 
elapsed time for this particular simulation is MAXTIME - TS 
= 3.333333E-5 seconds. 

The ZPRIME vector is created using values found in 
the TIME vector as shown in Equation 14. Each value of TIME 
is squared and multiplied by ce, Subtracting Z° from this 
product and taking the square root produces a ZPRIME ele- 
ment. Continuing with the previous example, two sample 
ZPRIME elements are 


ZPRIME(6) = SQRT(TIME(6)* * C% - Z%) = 1.304644E-2 


ZPRIME(7) = SQRT(TIME(7)* * C% - 27) = 1.848934E-2 
The algorithm loops until ZPRIME(STEP+1) through ZPRIME(64) 
have been calculated. } 

3. Generation of RHO2 array 

The next significant step is formation of the RHO2 

array. First the RHOl vector is created by calculating an 
increment and base value using the code 

INC = REAL(RHO) /REAL(NCOLS/2-1) 

BASE = REAL(-RHO1) - INC 
Assuming a value of 200 for RHO, INC = 6.451613 and BASE = 
=206.451613. This results in a RHO1l vector starting with 


-206.451613 in RHO1(1), with each successive element incre- 


mented by INC. This method ensures a value of 0.0 in 
RHO1(33), -200.0 (-RHO) in RHO1(2), and +200.0 (+RHO) in 
RHO1(64). As with the TIME vector, simply setting the RHOl 


vector to range from -RHO in RHO1(1) to +RHO in RHO1(64) 
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would result in identical values in RHO1(32) and RHO1(33). 
This would not meet the requirement that RHO1(33) contain 
singular information. 

After completing the RHO1l vector, its values are 
used to generate the RHO2 array. First the RHO1 vector 
values are placed in row 33 of RHO2 and its transpose in 
eolunne, 337 The program then fills in the blank array 
elements in quadrants II and III using the Pythagorean 
theorem with the values found in the respective row 33 and 
column 33 positions. Continuing with the example using RHO 
= 200, the value -206.451613 would be found in elements 
RHO2 (33,1) and RHO2(1,33) (column major format). Using the 
Pythagorean theorem, element RHO2(1,1) would receive the 
value 291.966671. This represents the value of RHO as 
calculated radially from the central position Or 
RHO2 (33,33). 

Since the RHO2 array is symmetrical, only quadrants 
It and III (left side) must be calculated. Although the 
time saved here is small, the algorithm using these values 
later in the program needs only the values found in the left 
half. Therefore the right side values are not calculated. 

4. Spatial excitation generation 

Generating the spatial excitation array is the next 
step. The first five selections represent single piston 
Spatial excitations centered about position (33,33) on a 64 


x 64 base array called PULSE. The five shapes available are 
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the square piston, circle piston, gaussian, pyramid, and 
truncated-pyramid. All five shapes require the user to 
specify the width (or diameter) of the excitation. The 
value-of width must be an odd integer to ensure symmetry 
about PULSE(33,33), and less than or equal to 63 since the 
base array size is 64 x 64. Small values between 9 and 15 
provide excellent results. | 

The square and circle selections create a simple 
spatial excitation with amplitude = 1.0. A resolution 
problem exists with the circular excitation. Simply, it is 
impossible to represent a circular shape using a_ small 
number of square elements. The larger the diameter, the 
smoother the outer surface becomes,’ but with a 64 x 64 array 
Size the range of available diameters is limited. The only 
solution is an increase in array size to 128 x 128 or 
larger. This is an excellent area for future expansion as 
micro-computers become faster and less memory limited. 

The Gaussian selection creates a circular-shaped 
excitation with amplitude following a gaussian distribution. 
The spatial excitation is truncated at the user specified 
diameter, and its width measured at the l1/e point. This 
excitation, as well as the pyramid and truncated-pyramid, 
represent non-piston sources with spatial variations in 
amplitude. 

The pyramid creates a four-sided tapered spatial 


excitation with a peak amplitude of 1.0 at the center and 
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tapering off until it equals 0.0 at the base. The size of 
the base is determined by the input value of width. Fous 
example, a pyramid excitation of width 11 would have its 
center at PULSE(33,33) and occupy the region bordered by 
rows 29 through 37 and columns 29 through 37. The amplitude 
would decrease linearly and equal zero at the outer edges of 
the above region. The truncated=-pyramid is similar except 
that the amplitude tapers off at a slower rate until it 
equals 0.0 at the array's outer edges. At the specified 
value of width, the pyramid is truncated, resulting in a 
"Square house" shape with tapered roof and square sides. 

The remaining selections all represent linear arrays 
of five elements each. An-individual element size of 5 x 5 
located on 10 unit wide center points was selected for 
convenience. Shape selections Lore the individual 
excitations are square, circular, and gaussian. Two 
additional arrays also use 5 square elements, but the 


amplitude of the elements are stepped or parabolic in shape. 


The stepped array has a center element of amplitude = 1.0, 
two adjacent elements of amplitude = 2/3, and two outer 
elements of amplitude = 1/3. The five elements thus create 
a stepped appearance: 1/3, 2/34 1.0j@)2/ 33 ee The 


parabolic array is similar to the stepped array, except the 
amplitudes follow a parabolic shape. 
After the spatial excitation is generated, the two 


dimensional Fast Fourier Transform is taken using the IMSL 
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routines FFT3D and FFTCC. The routine FFT3D computes the 
FFT of a complex-valued two-dimensional array by passing 
first the columns and then the rows as individual vectors to 
the routine FFTCC which computes the FFT of the vector argu- 
ment. The array returned by FFT3D is complex and has its 
peak information located at the outer four corners of the 
array, with the largest value in PULSE(1,1). The array is 
shifted to bring this peak information to the center using 
the subroutine SHIFT. This subroutine Swaps quadrant I with 
quadrant III and quadrant II with quadrant IV, thus trans- 
ferring the largest value from element PULSE(1,1) to element 
PULSE(33,33). This is why it is so important to arrange all 
arrays so their central value is in element (33,33). 
5. Filter Gp} for lossless media 

From this point the program branches based on the 
filter type selected. The algorithm for the lossless case, 
using filter Spi: is to take a value of ZPRIME, multiply the 
RHO2 array by this value and pass each product as the 
argument to the IMSL subroutine MMBSJN, where MMBSJN 
calculates the Bessel function of the first kind of order 
zero for a real number argument. The result returned from 
the Bessel function is multiplied by the respective element 
in array PULSE and stored in array WORK. The next value of 
ZPRIME is selected and the calculations again performed. 

As this portion of the program is the most time 


consuming, two specific actions were taken to optimize the 
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code. As previously explained, the number of rows 
identified by STEP are zeroed in the final array. Because 
of this, calculations start at row (STEP+1), thus saving the 
calculation of 4 rows for each iteration. Also, since the 
arrays RHO2 and PULSE are symmetrical in all four quadrants 
about element (33,33), only quadrant II needs to be 
calculated. With the default 64 x 64 array size and STEP = 
4, only 59 (64-STEP-1) subarrays of size 33 x 33 must be 
calculated for a total of 64,251 array elements. This 
compares favorably with the 262,144 array elements that 
would require calculation if the entire 64 x 64 array was 
calculated 64 times. It is important to note that the 
single most time-consuming process is the calculation of the 
Bessel function. 

After each 33 x 33 subarray is calculated, inline 
code is used to fill the entire WORK array. Quadrant II is 
flipped to fill quadrant III, then the entire left side is 
flipped to fill the right side. 

At this point the inverse two-dimensional Fast 
Fourier Transform of the WORK array needs to be taken. But 
instead of passing the entire array as an argument to a two 
dimensional inverse FFT subroutine, the array iS processed 
one column at a time. Using the IMSL subroutine FFT2C, 
which computes the inverse FFT (or FFT) of a complex-valued 
vector, each of the 64 columns are passed individually as a 


vector. After all the columns are transformed, only one 


28 


row, row 33, is passed as an argument. Assuming a symmetri- 
cal transducer, row 33 contains information from the central 
line across one axis of the observation plane. The inverse 
transform of this row represents the field values along this 
line for the specific value of ZPRIME, and it is this 
information that is used to sequentially build the final 
array called RESULT. The 59 iterations completed by the 
above algorithm will each supply one row, corresponding to 
an instant of time, in array RESULT. Specifically, with the 
default value of STEP, rows 5 (STEP+1) through 64 (NROWS) 
are filled. 

After completing row 64, the number of rows 
identified by STEP are set to zero as they do not contain 
any information. The entire RESULT array is then saved to 
disk, and the program terminates normally. 

6. Filter po for lossy media 

In the second case, using filter Sp2 for media with 
linear loss coefficient, the algorithm has two additional 
steps. The program implements Equation 17 and therefore 
must test each argument for RHO greater than AC/2. For 
RHO > AC/2, the argument is passed to the IMSL Bessel 
function subroutine MMBSJN as in the lossless case. But, if 
RHO < AC/2, the argument is instead passed to the IMSL 
subroutine MMBSIN which calculates the modified Bessel 
mimecion of ‘the first kind of order Zero. The result 


returned from either Bessel function is multiplied by the 
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respective element of PULSE. Then, as shown in Equation 17, 
the product is multiplied by an exponential term prior to 
storing in array WORK. This is the second difference in the 
lossy case. 

The algorithm loops, as in the lossless case, to 
fill the entire RESULT array. After saving the RESULT array 


on disk, the program terminates normally. 
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IV. NUMERICAL SIMULATION 


After careful debugging and testing on a point by point 
basis, the program was used to calculate the spatial impulse 
response for various source excitations. For each excita- 
tion shape, both the lossless and lossy case were calculated 
and the results plotted. The plots show the amplitude of 
the wave plotted against cross-direction and time at an ob- 
servation plane located above the X,Y plane. In all calcu- 
lations the height of the observation plane, Z, is 10 cm 
above the source plane. The value of A used for the lossy 
cases was 8.0E-3: This value was selected as a compromise 
and arrived at using trial and error. Smaller values of A 
caused only minor variations in shape while larger values 
caused gross distortion. All problems, unless noted other- 
wise, are calculated for a MAXTIME of 1.5E-4 seconds. Thus, 
with a Z value of 10 cm, each result runs from a starting 
time = Z/C = 6.666667E-5 seconds to 1.5E-4 seconds, a dura- 
tion of 8.333333E-5 seconds. 

As previously explained, all plots were obtained using 
the MATLAB MESH command. Two limitations exist with the 
MESH routine that require further explanation. The first 
limitation is,the way MATLAB sizes the plot. MATLAB uses 
input array values to fill a window of pre-determined size. 


This fixes the height of the plot so an excitation with 
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amplitude equal to 1.0 will plot identically) )co an 


excitation with amplitude other than 1.0. The second 
limitation builds on the first: MATLAB does not’ show 
numerically scaled axes when plotting an array of data. If 


it did, the different amplitudes for identically shaped 
plots would be apparent. 

An alternative way to obtain numerical information was 
found. Using the PLOT command, MATLAB has the ability to 
plot a vector in the X,Y domain with numerically scaled 
axes. Putting an array of data into PLOT creates a side 
view, plotting successive. "slices" of the array with time 
increasing in the positive X direction. Using this tech- 
nique, amplitude values could be obtained for a given plot. 
Since shape is the primary interest, amplitude information 
obtained in this manner was adequate. Three plotted arrays 
(Figures 9, 11, and 15) are included for illustration. 

Execution times were all less than the stated goal of 
thirty minutes, with an average of twenty-five minutes. 
These times were obtained using an AT&T 6300 micro-computer 
operating at 8 MHz. An 8087 numerical co-processor chip was 
installed and used to increase performance. Similar 
performance can be obtained with fast IBM XT compatibles or 
AT class micro-computers. The next generation of 
micro-computers utilizing the Intel 80386 microprocessor are 


the logical choice for continued program development. 
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A. CLASSICAL EXAMPLES 


The first two source shapes presented are the square 


piston and circular piston. 


eee 





Figure 7. Square piston spatial excitation. 


Figure 7 shows graphically the source excitation shape 
for a square transducer of width 11 and amplitude 1.0. The 
excitation is centered over the base array's central ele- 
ment, EBWESE (28yS3i)c In Figure 8 the calculated spatial im- 
pulse response for the square transducer in lossless media 


is shown, as is the array's orientation with time and space. 
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This orientation of time and space applies to all spatial 
impulse response plots. The origin of the time axis starts 
at Z/C where the potential is known to be a replica of the 
excitation source shape for lossless media. As time 
progresses, the potential is reduced until it forms two out- 


ward traveling "tails". As expected, the "tails" are square 
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Figure 8. Field for square transducer (lossless case). 


to match the spatial excitation shape, and slowly attenuate 
with time in the lossless case. Numerical information can 


be obtained from Figure 9 which is the side view of the 
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plot in Figure 8. The step function is clearly visible as 
is the jump in potential to 1.0 in row 5. From this point 
the potential peaks at a value slightly greater than 1.0 and 


then drops rapidly. At infinity the tails should reduce to 


zero. 
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Figure 9. Side view for lossless case. 


In Figure ‘10 the same excitation has been analyzed for 
the lossy case. Comparing the shapes in Figures 8 and 10 
reveals two distinct differences. The "tails" are no longer 


square but have taken on a triangular appearance. Larger 


aD 


values of A will produce a sharper triangular shape. The 
second major difference concerns the region between the 
tails. In the lossless case this region was at zero poten- 
tial. In the lossy case, the region does not approach zero 


but instead fills in, developing a shape of its own. 





Figure 10. Field for square transducer (lossy case). 


Another major difference is discovered when comparing 
Figures 9 and 11. Where the initial amplitude for the loss- 
less case was 1.0, the amplitude of the spatial excitation, 


the initial amplitude for the lossy case has been reduced to 
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approximately 0.30. 


This reduction in amplitude was ob- 


served for all test cases with a linear loss coefficient. 


The results obtained using a circular transducer are 


similar to those obtained with the square transducer. The 


spatial excitation for a circular transducer of width 15 and 


amplitude 1.0 is shown in Figure 12. Here a width of 15 was 
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Figure 11. Side view for lossy case. 


used to improve the resolution of the circular excitation. 


The spatial impulse response using this source is shown in 


Pigmre 13. 


Comparison of Figures 8 and 13 illustrates the 
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differences between the spatial excitation shapes. The ma- 
jor difference is the tails, which now have a circular 
cross-sectional shape. Some minor variations can also be 
observed in the large initial response early in time. 
Figure 14 shows the spatial impulse response, using the 
circular excitation of Figure 12, for the lossy case. The 


previous comments for the square transducer concerning the 





Figure 12. Circular piston spatial excitation. 


difference in amplitudes and shape of the "tails" for the 


lossless and lossy cases also apply to the circular 
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transducer. Figure 15 is a side view plot of Figure 14 
showing the reduction in amplitude is also present with the 
circular transducer. 

The results obtained when using the square-piston and 
circular-piston spatial excitation for both lossless and 
lossy media compare favorably with those obtained by Guyomar 
and Powers (Refs. ae Based on these results, it was 
determined that the program correctly implemented the 
propagation models, and therefore could be used to compute 


the spatial impulse response uSing new types of excitation. 
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Figure 13. Field for circular transducer (lossless case). 
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B. NEW EXCITATION SHAPES 

The program has to ability to analyze any excitation 
source shape that will fit on the 64 x 64 base array. Two 
new single excitation shapes and four multiple element exci- 


tation shapes are in the program. The results for three of 





Figure 14. Field for circular transducer (lossy case). 


these, the pyramid, the 5 element linear array with constant 
amplitude, and the 5 element linear array with stepped am- 


plitude are included. 


A pyramid shaped spatial excitation of width 11 and am- 
plitude 1.0 is shown in Figure 16. The computed spatial im- 
pulse response for this excitation in lossless media is 


shown in Figure 17. Since this particular waveshape has a 
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Figure 15. Side view for lossy case. 


low spatial frequency content, the shape of the wave stays 
much the same for small values of time. Figure 18 shows the 


spatial impulse response for the same pyramid excitation in 


lossy media. 
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Figure 19 shows the 5-element linear array configura- 


tion. For this array the elements are each 5 x 5 with their 


centers positioned 10 units apart. This produces a space of 


5 units between each element. The calculated spatial im- 


pulse response for this excitation is shown in Figure 20. 





Figure 16. Pyramid shaped spatial excitation. 


The shape for each individual excitation is similar to the 
Single excitation response in Figure 8 until the tails start 


to spread out and interfere with each other. In the lossy 
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case, Figure 21, the pattern is similar but with increased 


interference in the region between the individual responses. 





Figure 17. Field for pyramid excitation (lossless case). 


The spatial excitation for the last test case presented 
is shown in Figure 22. This is a 5-element linear array 
Similar to Figure 19 but with stepped amplitude. The center 
element has amplitude of 1.0 as before. The two elements 
adjacent to the center have amplitudes of 2/3, and the outer 
two elements have amplitudes of 1/3. This particular case 


is presented to show the flexibility of the program. Figure 
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23 shows the spatial impulse response in lossless media. 
The results are Similar to those of Figure 20 except for the 
dominance of the center element. Figure 24 shows the spa- 
tial impulse response for the same excitation in lossy me- 


dia. 





Figure 18. Field for pyramid excitation (lossy case). 
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Field for five element linear array 


(lossless case). 


Figure 20. 
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Figure 21. Field for five element linear array 
(lossy case). 
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Figure 22. Excitation for five element linear array 
with stepped amplitude. 
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Figure 23. Field of five element linear array with 
stepped amplitude (lossless case). 
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Field for five element linear array with 


stepped amplitude (lossy case). 


Figure 24. 
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V. SUMMARY 


This thesis has investigated the ability to perform 
complex simulations of scalar wave propagation through loss- 
less media and media with a loss coefficient linear in fre- 
quency on a micro-computer. A program was written in 
Microsoft Fortran V3.31 which allows the user to select or 
specify the initial: conditions and media type, then calcu- 
lates the resulting spatial impulse response. The set goal 
of thirty minutes for each simulation was achieved through 
code optimization and partial array calculation permitted by 
symmetry. 

After initial program verification, the spatial impulse 
responses for both square-piston and circular-piston spatial 
excitation were computed. The results obtained agreed with 
previously accepted results. The spatial impulse response 
was then calculated for three new spatial excitation 
sources, a pyramid-shaped excitation, a five element linear 
array excitation with constant amplitude, and a five element 
linear array excitation with stepped amplitude. The 
computed spatial impulse response for all simulations have 
been plotted and are presented for visual interpretation. 

Future development should concentrate in two areas. 
First, the third model for media with a loss coefficient 


which is quadratic in frequency must be added to the 


or 


program. This will provide all the tools necessary for 
analysis of wave propagation through air, liquid, gas, and 
biological tissue. Second, the base array size must be in- 
creased beyond the present 64 x 64. An array size of 64 x 
64 proved adequate for program development and testing but 
increased size, and therefore increased resolution, is nec- 
essary for serious analysis and evaluation. The limitation 
in all presently available micro-computers which restricts 
array size is processing speed. For example, increasing ar 
ray size to 128 x 128 would increase execution time by a 
factor of four, requiring almost two hours for each simula- 
C2LOn: run. Transferring program development to one of the 
newer micro-computers utilizing the Intel 80386 micro- 
processor will allow larger array sizes while maintaining 


reasonable execution times. 
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APPENDIX: FORTRAN SOURCE CODE 





This appendix contains the FORTRAN source code for the 


program developed in the thesis. This program was compiled 


and run using Microsoft Fortran Version 3.31. 


SLARGE i: 
SNOFLOATCALLS 


oO a a ae Ono OO A aaa I OOo oe fe Aaa a a ea ea A ea eae ae a 


Here te te tek te te te tected dete tote te deve tek terete tke tick keiknekkKKKKKKKKKKKKKKKKKKKKKRK KKK KK RK 


COMMENTS : 
THIS PROGRAM WAS WRITTEN BY LT. TIMOTHY MERRILL, USN 


IT IS WRITTEN IN MICROSOFT FORTRAN V3.31 AND IS CAPABLE OF RUNNING 
ON ANY IBM OR COMPATIBLE MICRO-COMPUTER WITH LITTLE OR NO 
MODIFICATION. THE USE OF A 8087 NUMERICAL CO-PROCESSOR IS HIGHLY 
RECOMMENDED TO REDUCE EXECUTION TIME. 


TO REDUCE EXECUTION TIME THE PROGRAM IS WRITTEN TO USE A 64 X 64 
BASE ARRAY SIZE, AND IS WRITTEN AS SINGLE PRECISISON. FUTURE 
MODIFICATION TO A 128 X 128 (OR LARGER) ARRAY SIZE AND/OR DOUBLE 
PRECISION IS VERY EASY. 


THE PROGRAM CALLS THE FOLLOWING SUBROUTINES WHOSE SOURCE CODE IS 
NOT INCLUDED HERE: 


FFT2C IMSL ROUTINE TO COMPUTE FFT OF A COMPLEX 
VALUED SEQUENCE OF LENGTH EQUAL TO POWER 
OF TWO 
FFT3D IMSL ROUTINE TO COMPUTE THE FFT OF A COMPLEX 
VALUED 1,2, OR 3 DIMENSIONAL ARRAY 
FFTCC IMSL ROUTINE CALLED BY FFT3D 
MMBSJIN IMSL ROUTINE TO CALCULATE BESSEL FUNCTION OF 
THE FIRST KIND OF ZERO ORDER WITH REAL ARGUMENT 
MMBSIN IMSL ROUTINE TO CALCULATE MODIFIED BESSEL 


FUNCTION OF THE FIRST KIND OF ZERO ORDER WITH 
REAL ARGUMENT 


T.D.M. 20 JANUARY 1987 


We ve Ye Hee ve te ke ve te te ve te te Ke te ve Ye te te te te te te te ke ke He Ve He Ve Ve Ve Ve Ve Ye Ve He Ve Fe ve He Xe Ve Ve Ve Ve Ye Ve Fe ye Ve Ve Fe Te ve Ye de re Ve re Ve Ve ve ve He 


REAL A, AZ, Bye, BASE, BI,BR,C,CZ,CENTER 
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Oo n 2 oan aA a Ona a 040 a Oa AO a 4 


REAL 
REAL 


COMPLEX 
COMPLEX 


INTEGER 
INTEGER 
INTEGER 


DIMENSION 
DIMENSION 
DIMENSION 


CHARACTER 
CHARACTER 
CHARACTER 


HWIDTH, INC ,MAXTIME , MEANX, PI 

R123 , RESULT ,RHO1, RHO2, RWK, SIGMA, SIGMAZ 
TEMP1,TEMP2 , TEMP3 , TEMP4 , TEMPS 

TIME ,TS,IWOPI ,X,Y,Z, 22, 2PRIME, TVECIA 


C123 ,CWK,CTEMP1 ,CTEMP2 ,CTEMP3 
PULSE , WORK 


COL, FTYPE,H,1I, IER, IJOB, IWK, IWK1, J 
M,N,NAMLEN ,NCOLS, NROWS ,R1,R2 
ROW, RHO, SHAPE , STEP ,WIDTH, XPOS, YPOS 


CWK(64) , IWK( 534) , IWK1(7) ,RHO1(64) , RWK(534) 
PULSE (64,64), RESULT (64,64) ,RHO2(64, 64) 
TIME (64) ,WORK(64,64) , ZPRIME(64) 


FN*S, FNAME*3, 11,12 
EXT1*4 , EXT2*4, EXT3*4, EXT4*4 
FNAME1*12, FNAME2*12, FNAME3*12, FNAME4*12 


We te Ye He Ye He He He ve Fe He He He He Ye He He we He Ye He ve He He He He Ye He He He He He He He He He He He He ve He Ve He He ve We Ye He He He Ye Fr Ve He ve He He He He He He He ve He He 


CONSTANTS 


We ve ve ve ve ve He He ve tr te ve ve ve ve ve ve ve Ye ve We ve ve We ve ve ve te ve We ve ve We ve ve ve We We We ve He We ve ve ve ve He ve We ve ve ve te ve ve ve ve ve We ve or ve ve tr ve 


C 
C2 
H 
M 


MEANX 
N 
NCOLS 
NROWS 
SIGMA 
SIGMA2 
STEP 
TWOPI 


SPEED OF SOUND IN VACUUM (1500 METERS / SEC} 

SPEED OF SOUND SQUARED (METERS**2/SEC**2) 

HALF OF NROWS (OR NCOLS) PLUS 1 

INTEGER USED TO SPECIFY SIZE OF VECTOR SENT TO 
FFT2C (REFER TO IMSL SOURCE CODE FOR FFT2C) 

USED FOR GAUSSIAN IMPULSE FUNCTION 

USED ON BESSEL FUNCTION CALLS TO SIGNIFY ZERO ORDER 

DEFAULT NUMBER OF COLUMNS 

DEFAULT NUMBER OF ROWS 

USED FOR GAUSSIAN IMPULSE FUNCTION 

SIGMA SQUARED 

NUMBER OF ROWS SET TO ZERO - SIMULATES STEP FUNCTION 

2: tet 


We de tre He He ve He de ve He ve te te Fe ve ve ve We We ve He re Ve We He te He We Hc We Ve ve He ve We He Fe We ve He We He He He ve He He He We He ve He He We ve He He ve He He He He He He 


C 
Nae 
M=6 
STEP = 4 
NROWS 
NCOLS 


i 


1500.0 


64 
64 


PI. = 3.14159265 


SIGMA = 1.0 
MEANX = 0.0 
CZe= Cle iC 


H = NROWS/2 + 1 
TWOPL, #* 2.0 * PI 
SIGMA2 = SIGMA * SIGMA 
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aQ 


He He He He HAH HH HHH I HTK IH HK HTK III HK EKER EKER HERE HH KK 


GET USER INPUT 

Fe Fe HK KH eH KK KKH HI HK IKKE KEK EKER KEKE KHAKI KKK KKK 
XPOS = 1 

YPOS = 2 

CALL CLRSCR 

CALL GOTOXY(1, 1) 

WRITE(*,*) ‘ENTER REQUESTED VALUES. RECOMMENDED VALUES IN ( )? 
CALL GOTOXY (YPOS , XPOS) 

WRITE(*,*) ‘ENTER FILENAME (8 CHAR MAX): ’ 
CALL GOTOXY (YPOS , XPOS+32) 

READ(*,’(A)’) FN 

CALL GOTOXY (YPOS+1,XPOS) 

WRITE(*,*) ’SELECT FILTER TYPE’ 

CALL GOTOXY (YPOS+2 ,XPOS+8) 

WRITE(*,*) ’<1> GP1 (LOSSLESS MEDIA)? 
CALL GOTOXY (YPOS+3, XPOS+8) 

WRITE(*,*) ’<2> GP2 (LOSSY MEDIA)’ 

CALL GOTOXY (YPOS+4 , XPOS) 

WRITE(*,9) ’ ENTER 10R 2: ’ 

READ(*,*) FTYPE 

CALL GOTOXY (YPOS+5 , XPOS) 

WRITE(*,*) ’SELECT SOURCE SHAPE’ 


“CALL GOTOXY (YPOS+6 , XPOS+4) 


WRITE(*,*) ’<1> SQUARE’ ©~ 
CALL GOTOXY (YPOS+7 ,XPOS+4) 
WRITE(*,*) ’<2> CIRCLE’ 
CALL GOTOXY (YPOS+8 , XPOS+4) 
WRITE(*,*) ’<3> GAUSSIAN’ 


- CALL GOTOXY (YPOS+9 ,XPOS+4 ) 


WRITE(*,*) *<4> PYRAMID’ 

CALL GOTOXY (YPOS+10,XPOS+4) 

WRITE(*,*) *<5> TRUNCATED PYRAMID’ 

CALL GOTOXY (YPOS+11,XPOS+4) 

WRITE(*,*) ’<6> 5 PULSE LINEAR ARRAY (SQ.)’ 
CALL GOTOXY (YPOS+12 , XPOS+4) 

WRITE(*,*) *’<7> 5 PULSE LINEAR ARRAY (GAUSS. )’ 
CALL GOTOXY (YPOS+6 , XPOS+40) 

WRITE(*,*) *’ <8> 5 PULSE LINEAR ARRAY (STEPPED)’ 
CALL GOTOXY(YPOS+7 , XPOS+40) 

WRITE(*,*) *’ <9> 5 PULSE LINEAR ARRAY (PARABOLA)’ 
CALL GOTOXY (YPOS+8 , XPOS+40) 

WRITE(*,*) ’<10> RESERVED FOR FUTURE USE’ 

CALL GOTOXY(YPOS+tS9 , XPOS+40) 

WRITE(*,*) ’<11> RESERVED FOR FUTURE USE’ 
CALL GOTOXY (YPOS+10 , XPOS+40) 

WRITE(*,*) ’<12> RESERVED FOR FUTURE USE’ 
CALL GOTOXY (YPOS+11,XPOS+40) 

WRITE(*,*) ’<13> RESERVED FOR FUTURE USE’ 
CALL GOTOXY(YPOS+12 ,XPOS+40) 

WRITE(*,*) *<14> RESERVED FOR FUTURE USE’ 
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CALL GOTOXY (YPOS+13 , XPOS ) 

WRITE(*,9) ’ ENTER SELECTION: ’ 

READ(*,*) SHAPE 

IF (SHAPE .LT. 6) THEN 
CALL GOTOXY (YPOS+14 , XPOS) 
WRITE(*,9) ’ ENTER SOURCE SIZE (ODD INTEGER <= 63) (11): 
READ(*,*) WIDTH 

ENDIF 

CALL GOTOXY(YPOS+15, XPOS) 

WRITE(*, *)' “ENTER RHOG™3(2Z00>: 

CALL GOTOXY(YPOS+15 , XPOS+26) 

READ(*,*) RHO 

CALL GOTOXY (YPOS+16 , XPOS ) 

WRITE(*,*) “ENTER Z (METERS) (0.1)3= 

CALL GOTOXY (YPOS+16 , XPOS+26 ) 

READ(*,*) Z 

CALL GOTOXY (YPOS+17 , XPOS ) 

MERITE(*,7) 2/C 

CALL GOTOXY(YPOS+17, XPOS+44 ) 

READ(*,*) MAXTIME 

IF (PL YPE EQ. 2) TEEN 
CALL GOTOXY (YPOS+18, XPOS ) 
WRITE(*,*) ’ENTER VALUE FOR A (0.008): ’ 
CALL GOTOXY (YPOS+18 , XPOS+265 ) 


READ(*>*)% 
AZ = A * A 
ENDIF 


FORMAT(’ ENTER MAXTIME (REAL > ’,1PE14.6,’):’) 
FORMAT (A\ ) 

HWIDTH = REAL(WIDTH)/2 

Z2=ezZ*Z 

TS = Z/C 

CALL CLRSCR 

CALL GOTOXY(1,1) 

WRITEC*,*) *CALCULATING..: <” 


, 


Fete ve te te We ve ve He de He ve He ve ve He ve ve We ve ve He te ve He ve ve He He He ve ve ve He te ve He ve He ve ve te ve He ve He He ve He He ve He ve ve He He ve He ve ve ve ve He He He 


INITIALIZE ARRAYS AND VECTORS TO ZERO 


He te He ve te He He ve te ve ve He ve ve ve te ve ve ve ve ve ve ve ve ve He te He ve ve ve He ve He He Ve Ye Ye we ve ve ve ve He ve He ve He He He He ve ve We ve We Ye ve ve ve ve ve He ve He 


DO 5 ROW = 1,NROWS 
TIME (ROW ) 0.0 
RHO1 (ROW) 0.0 
DO 5 COL = 1,NCOLS 
PULSE(COL,ROW) = (0.0,0.0) 
WORK(COL,ROW) = (0.0,0.0) 
RESULT(COL,ROW) = 0.0 
RHOZ(COL,ROW) = 0.0 
CONTINUE 
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qgqaagaaaaaaaaaa a 


10 


qa aaa 


He He Fe te HHH HICH HH HK TATA IITHTHTITATATIAHE KITTTKATEATHIREEK I IIK 


OPEN OUTPUT FILES 


He He HK HH HHH HHH HHI IKI KIKI KKK KEKE REE KITE HIKE KEKKK KK 


PROGRAM GETS NAME OF PROBLEM <FN> FROM USER AND FORMS: 


FILENAME. EXT FILE DESCRIPTION 
<FN>. IMP ARRAY OF INPUT IMPULSE FUNCTION 
<FN>. TXT SUMMARY INFORMATION FOR THIS PROBLEM 
<FN>.DAT FINAL RESULTS ARRAY 


NOTE 1: ALL FILES ARE IN ASCII TEXT FORMAT 

Ye ve He Ye te He Ye He Fe te He i te He te He ee He He Fete We He He He He He He He Pete Fe te Ye He He Fete He He He We He Fe re Fe He He He ve We ie He Fe He eH He Weve He 
EXT1 = ’. IMP’ 

EXT2 = ’ TXT’ 

EXT3 = ’.DAT’ 

NAMLEN = 0 

Ts’ ° 

DO 10 I=1,8 

T2 = FN(I:I) 

IF (T2 .NE. T1) THEN 
FNAME(I:1) = FN(I:1) 
NAMLEN = NAMLEN + 1 

ENDIF 

CONTINUE 

FNAME1(1:NAMLEN) = FNAME 
FNAME 1 (NAMLEN+1 : NAMLEN+4 ) 
FNAME2(1:NAMLEN) = FNAME 
FNAME2 (NAMLEN+1 : NAMLEN+4 ) 
FNAME3(1:NAMLEN) = FNAME 
FNAME3 (NAMLEN+1 : NAMLEN+4 ) 
IF (NAMLEN .LT. 8) THEN 

FNAME1(NAMLEN+5:12) => 

FNAME2 (NAMLEN+5 : 12) 

FNAME3 (NAMLEN+5 : 12) 

ENDIF 

OPEN(1,FILE=FNAME1, STATUS=’ NEW’ ) 
OPEN(2, F ILE=FNAME2 , STATUS=' NEW’ ) 
OPEN (3, FILE=FNAME3 , STATUS=’ NEW’ ) 


EXT1 


EXT2 


EXT3 


3 3 


3 r 


He ve He He Ye Fe Ye He Fe Fe He He He Fe Ye He Ye He He He He He He Ye Ye ve He He He te He Fe He Ye Ye Ye Fe Ye Ve He He We Ye Fe He He He He He He We We He He He He He Ye Ye He He He He He 
WRITE HEADER INFORMATION TO <FN>.TXT 
We tote He te ie i ve He te He te He Fe He Fe He He He We He He te ie ve He Ye Fe Fe Fe ie te ve Ye We He He He Ye He Ye He Ye Fe He Ye Fe He He He He He He He ve He He ve Ye Fe ve Fe He Yee 
WRITE(2,11) FN 
ite CPIYPE .EQ. 1) THEN 
WRITE(2,*) *® FILTER IS GP1 FOR LOSSLESS MEDIA’ 
ELSELF (FIYPE .EQ. 2) THEN 
WRITE(Z,*) ’ FILTER IS GP2 FOR LOSSY (LINEAR) MEDIA’ 
ENDIF 
TE CSHAPE .EQ. 1) THEN 
WRITE(2,12) WIDTH,NCOLS , NROWS 
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el 
12 


13 


14 


114 


LS 


i 
16 
Le 
18 
19 


ELSEIF (SHAPE .EQ. 2) THEN 


WRITE(2,13) WIDTH,NCOLS , NROWS 


ELSEIF (SHAPE .EQ. 3) THEN 


WRITE(2,14) WIDTH, NCOLS , NROWS 


ELSEIF (SHAPE .EQ. 4) THEN 


WRITE(2,111) WIDTH,NCOLS,NROWS 


ELSEIE (SHAPE (20275) 1cEN 


WRITE(2,112) WIDTH,NCOLS,NROWS 


ELSEIF (SHAPE .EQ. 6) THEN 
WRITE(2,113) 

ELSEIF (SHAPE .EQ. 7) THEN 
WRITE(2,114) 

ELSEIF (SHAPE .EQ. 8) THEN 
WRITE(2,115) 

ELSEIF (SHAPE .EQ. 9) THEN 
WRITE(2,116) 

ENDIF 

WRITE(2,15) STEP 

WRITE(2,16) 2 

WRITE(2,17) RHO 

WRITE(2,18) MAXTIME 

IF (FTYPE .EQ. 2) THEN 
WRITE(Z,19) A 

ENDIF 


FORMAT (’ SUMMARY INFORMATION FOR ’ ,A/) : 

FORMAT (’ IMPULSE FUNCTION IS A SQUARE WITH WIDTH ’,I2,’ ON A BAS 
TE OR Lay a 4 oo 

FORMAT (’ IMPULSE FUNCTION IS A CIRCLE WITH DIAMETER ’,I2,’ ON A 


+BASE OF ,13,7  X*, 13) 
FORMAT (’ IMPULSE FUNCTION 
BASE OF’ ,13,° XX’ ,13) 


IS 


FORMAT(’ IMPULSE FUNCTION IS 
+E OF’,I3,’ X’,13) 
FORMAT(’ IMPULSE FUNCTION IS 


+ON A BASE OF 713)" XxX’, 13) 


GAUSSIAN WITH DIAMETER ’,1I2,’ ON A 
PYRAMID WITH WIDTH ’,I2,’ ON A BAS 
TRUNCATED PYRAMID WITH WIDTH ’,12,’ 
5S PULSE LINEAR ARRAY WITH SQUARE PU 
S PULSE LINEAR ARRAY WITH GAUSSIAN 
5S PULSE LINEAR ARRAY WITH STEPPED A 


5 PULSE LINEAR ARRAY WITH PARABOLIC 


’ SECONDS’ ) 


FORMAT (’ IMPULSE FUNCTION IS 
+LSES’ ) 

FORMAT (’ IMPULSE FUNCTION IS 
+PULSES* ) 

FORMAT (’ IMPULSE FUNCTION IS 
+MPLITUDE’ ) 

FORMAT ¢’ IMPULSE FUNCTION IS 
+ AMPLITUDE’) 

FORMAT ¢ ’ SIEPSOILZE = * , 12) 
FORMAT (’ Su=" PETS. 7," METERS” ) 
FORMAT(’ RHO = ’,13) 

FORMAT (’ MAXTIME = ’,1PE14.7, 
FORMAT ¢ ’ A = ’,1PE14.7) 
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a 


qaaana 


21 


22 


qgqaQagaaaaa 


30 


aaa a a a 


50 
C 


60 
C 


70 
C 


Fe He He He He HHH IH He eK KKK KEK KKKKAKAKHKHKEK AK KH eK ee tk ee 


INITIALIZE TIME VECTOR 
Fe HHI IOI HITCH ICT IITA IT IIIA IIT TI IT III I TA TTT TTI IIIA AAI III IK 
INC = (MAXTIME-TS) /REAL(NROWS-STEP-1) 
DO 20 I = (STEP+1),NROWS 
TIiMeECh ests + (I-STEP-1) * INC 
CONTINUE 


He We He ie He He He He He He He ie ee He He He He He He HH HHH KKK KH KKHKHKKKKKEKKKKKKKKENKKKKKKKKK KK 


OUTPUT TIME SUMMARY TO <FN>.TXT 

IACI IRI ISI III I IOI OIC IOI IOI RII IOI TIO ITI TON IT HR II 
WRITE(2,*) 

WRITE(2,9) ’TIME SUMMARY: ’ 

WRITE(2,21) TS,STEP+1 

WRITE(2,22) INC 

FORMAT (’ ZERO TIME STARTS AT T=Z/C =’,1PE14.7,’ SECONDS IN ROW(’ 
aie, )” ) 

FORMAT (’ AND INCREMENTS BY ’,1PE14.7,'’ SECONDS PER ROW’/) 


HECK K KKK KK KKK KKK KKK KKK KKEKKKKKKKKKKKKKEKKKKKKK 


CREATE Z-PRIME VECTOR: SQUARE EACH VALUE OF TIME, 
MULTIPLY BY C SQUARED AND SUBTRACT Z SQUARED. THEN 
TAKE SQUARE ROOT. 
He eH He Hee NE IO III ICH IOT ICT IEICE IE IIIT DIINO II IT IIIA IIIA ICT IIIT ICT Te HICK tok 
DO 30 I = (STEP+2),NROWS 
ZPRIME(I) = SQRT((TIME(I) * TIME(I) * C2) - 22) 
CONTINUE 


He He He He He He He He He He He He He Me He Ye He te He He He He Me He He He He He He He Ye ve He He He He He Ye He Ye He He He ve ve He We He He He He He He He He ee ke eke ke ke ke ee 


INITIALIZE RHO VECTOR AND GENERATE 33 X 64 RHO MATRIX 
- ONLY ONE HALF OF THE MATRIX HAS TO BE CALCULATED SINCE 
THE OTHER SIDE IS SYMMETRICAL 
We He He He He Ye He He He He He He He He He He He He He te Hee He te He te he ee Ke he He te He He He Hee He He He He He He He He He He He ie Ke te ete ee ee He He 
INC = REAL(RHO) /REAL(NCOLS/2-1) 
BASE = REAL(-RHO) - INC 
~RHO1(1) = BASE 
DO 50 I = 1,NCOLS-1 
RHO1(I+1) = BASE + I * INC 
CONTINUE 


DO 60 I = 1,NROWS 
RHO2(H,I) = RHO1(I) 
RHO2(1I,H) = RHO1(T) 

CONTINUE 


DO 70 ROW = 1,NROWS 
DO 70 COL = 1,H 
RHO2(COL,ROW) = SQRT(RHOZ(COL,H)**2 + RHO2(H,ROW)**2) 
CONTINUE 


ao 


G RKEKKKKKKKKKKKKKEKKKKKKEKKKKKEKEKKKEKKKKKKEKKKKKKEKKKKEKKKEKKKKKKKKKKK 


C OUTPUT RHO SUMMARY TO <FN>.TAT 
c eH Pe RFR IH ICT ICICI ICICI FEIT IE ICT IE ICT TOTO IT IIT AIT AE IAA IA AIA IA IAA AIA AIA AAA AIAN 
WRITE(2,9) ’RHO SUMMARY: ’ 
WRITE(2,44) H,NCOLS/2+1 
WRITE(2,41) BASE 
WRITE(2,42) INC 
WRITE(2,43) RHO2(1,1) 


41 FORMAT (’ BASE VALUE = (-RHO - INC) = ’ IPE14.7) 
42 FORMAT (’ INCREMENT = RHO/(64/2 - 1) eT AEE LS 
43 FORMAT (’ MAXIMUM VALUE IS AT RHO2(1,1) = ’,1PE14.7) 


44 FORMAT (’ ARRAY IS CENTERED ABOUT ("13,5 ,412, 9.) 


KEK KKKKKKKKKKKEKKKKKKKKKKKKAKEKKKKKKKKKRKKKKAKKEKKKKKKKKEKKKKKEKKKKK 


C 
C 
C GENERATE EXCITATION - CENTER POINT MUST BE AT (33,33) 
C 
C 
C 


RKAKKAMKAKKAKKKKKKKKKKKMKKKKKKKKKKKEKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK 


axe GENERATE SQUARE EXCITATION > 
IF (SHAPE .EQ. 1) THEN 
DO 1200 ROW = (NROWS/2+1-WIDTH/2), (NROWS/2+1+WIDTH/2) 
DO 1200 COL = (NCOLS/2+1-WIDTH/2) , (NCOLS/2+1+WIDTH/2) 
PULSE(COL ,ROW) = (1.0,0.0) 


1200 CONTINUE 
C 
c «xs GENERATE CIRCULAR PULSE *** 


ELSEIF (SHAPE .EQ. 2) THEN 
DO 1300 ROW = (NROWS/2+1-WIDTH/2), (NROWS /2+1+WIDTH/2) 
DO 1300 COL = (NCOLS/2+1-WIDTH/2), (NCOLS/2+1+WIDTH/2) 
IF (SQRT(REAL(COL-H)**2+REAL(ROW-H)**2) .LT. HWIDTH) THEN 
PULSE(COL,ROW) = (1.0,0.0) 


END IF 
1300 CONTINUE 
C 
C am GENERATE CIRCULAR EXCITATION WITH GAUSSIAN AMPLITUDE =~ 
ELSEIF (SHAPE .EQ. 3) THEN 
TEMP1 = 1/SQRT ( TWOPI*SIGMA2) 
TEMP2 = 1/(2*SIGMAZ) 
DO 1400 ROW = (NROWS/2+1-WIDTH/2), (NROWS/2+1+WIDTH/2) 
DO 1400 COL = (NCOLS/2+1-WIDTH/2) , (NCOLS/2+1+WIDTH/2) 
TEMP3 = SQRT(REAL(COL-H)**2+REAL (ROW-H) **2) 
IF (TEMPS .LI. HWIDTH) THEN 
PULSE(COL,ROW) = TEMP1*EXP(-TEMP2*(TEMP3-MEANX )**2) 
ENDIF 
1400 CONTINUE 
C 
C aay GENERATE PYRAMID EXCITATION ih 


ELSEIF (SHAPE .EQ. 4) THEN 
Rl = (NROWS/2) - (WIDTH/2) 
R2 = (NROWS/2 + 2) + (WIDTH/2) 
ING =.0/ (WIDTH 2-3 1) 
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DO 1500 I = 1,WIDTH/2 
DO 1500 ROW = R1+I,R2-I 
DO 1500 COL = R1+I,R2-I 
PULSE(COL,ROW) = INC * I 
1500 CONTINUE 
PULSE(H,H) = 1.0 


c *** GENERATE TRUNCATED PYRAMID EXCITATION *** 
ELSEIF (SHAPE .EQ. 5) THEN 
Rl = (NROWS/2) - (WIDTH/2) 
R2 = (NROWS/2 + 2) + (WIDTH/2) 
INC = 1.0/(NCOLS/2) 
BASE = (NCOLS/2 - WIDTH/2 + 1) * INC 
INC = 1.0/H 
DO 1600 I = 1,WIDTH/2 
DO 1600 ROW = R1+I,R2-I 
DO 1600 COL = R1+I,R2-I 
PULSE(COL,ROW) = BASE + INC * I 
1600 CONTINUE 
PULSE(H,H) = 1.0 


E 
Cc **x* GENERATE 5 ELEMENT LINEAR ARRAY - SQUARE *** 
ELSEIF (SHAPE .EQ. 6) THEN 
WIDTH = 5 
DO 1700 I = 1,WIDTH 
COL = I * 10 
DO 1700 J = 1,WIDTH 
COL = COL + 1 
DO 1700 ROW = (H - WIDTH/2),(H + WIDTH/2) 
PULSE(COL,ROW) = (1.0,0.0) 
1700 CONTINUE 
c 
C *#** GENERATE 5 ELEMENT LINEAR ARRAY - GAUSSIAN *** 
ELSEIF (SHAPE .EQ. 7) THEN 
WIDTH = 5 
TEMP1 = 1/SQRT(TWOPI*SIGMA2) 
TEMP2 = 1/(2*SIGMA2) 
DO 1800 ROW = (NROWS/2+1-WIDTH/2), (NROWS /2+1+WIDTH/2) 
DO 1800 COL = (NCOLS/2+1-WIDTH/2), (NCOLS/2+1+WIDTH/2) 
TEMP3 = SQRT(REAL(COL-H)**2+REAL (ROW-H) **2) 
IF (TEMP .LT. REAL(WIDTH/2)) THEN 
TEMP4 = TEMP1*EXP(-TEMP2* ( TEMP3-MEANX ) **2) 
PULSE(COL-20,ROW) = TEMP4 
PULSE(COL~-10,ROW) = TEMP4 
PULSE (COL , ROW) = TEMP4 
PULSE(COL+10,ROW) = TEMP4 
PULSE (COL+20,ROW) = TEMP4 
ENDIF 
1800 CONTINUE 
Cc 
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C eal GENERATE 5 ELEMENT LINEAR ARRAY - STEPPED a 
ELSEIF (SHAPE .EQ. 8) THEN 
WIDTH = 5 
TEMP1 = 1.0/3.0 
TEMPZ “2.0/3.0 
TEMP3 = 1.0 
DO 1900 COL = 1,WIDTH 
DO 1900 ROW = (H-WIDTH/2), (H+WIDTH/2) 
PULSE (COL+10 , ROW) TEMP1 
PULSE (COL+20 , ROW) TEMP2 
PULSE (COL+30 , ROW) TEMP3 
PULSE (COL+40 , ROW) TEMP2 
PULSE (COL+50 , ROW) TEMP1 
1900 CONTINUE 


wae LINEAR ARRAY OF 5 ELEMENTS - PARABOLIC AMPLITUDE ane 

IMPLEMENTS THE FOLLOWING FORMULA: 

(KX ="H)*"2 we -465"% Fo * (Y=) 
WITH 
H.=.0 
256 
K = 256 

TO GIVE A DOWNWARD OPENING PARABOLA THAT HAS A PEAK OF 1.0 
AND IS EQUAL TO ZERO AT X = -32, 32 (ROWS 2 AND 64). 
EQUATION IS EVALUATED AT X = 0, 10 AND 20. 


(1-0): 0-0 Oe 0 OD) Oe OO) a 
ro 
i} 


ELSEIF (SHAPE .EQ. 9) THEN 
“WIDTH = 5 
TEMPL. = (4.%250. = 20%"2)/ (4. "2562 
TEMP2 (4.*256. - 10**2)/(4.*256. ) 
TEMPSE= "04.%290:. = -0**%Z)/ (4.."256... 
DO 2000 COL = 1,WIDTH 
DO 2000 ROW = (H-WIDTH/2), (H+WIDTH/2) 


PULSE(COL+10,ROW) = TEMP1 
PULSE (COL+20,ROW) = TEMP2 
PULSE (COL+30,ROW) = TEMP3 
PULSE (COL+40,ROW) = TEMP2 
PULSE (COL+50,ROW) = TEMP1 
2000 CONTINUE 
C 
ENDIF 
S 
E ZEIGE OICICICIIIOIICIOIOIOIICISIOI ICICI ICID IOI GOI ICICI IOI IOIOI IOI ICAI ISTO IOI I Tt tote tri 
C OUTPUT IMPULSE FUNCTION ARRAY TO <FN>. IMP 
C TOIAIRIIICIOICI ICICI IOI CIEIFIOI ICISIOISICICIOI IOI GOIOIIIOIIOII IOI ICIS IOI II III IOI ICI io 


DO 3000 I = 1,NROWS 
WRITE(1,102) (REAL(PULSE(J,1I)),J = 1,NCOLS) 
3000 CONTINUE 
CLOSE(1) 


62 


a a 


qAaaaa 


gQaangnaa a a 


210 


qAaana 


220 


qaqgaaa 


215 


MHA KEKEKEERRREREEREKKKKKKKKKKAKKRAK 


PASS PULSE TO 2-D FFT SUBROUTINE. IJOB = 1 MEANS TAKE FFT 
FAI II III II IT II TO IT ITI TIT IITA IIIA TAI AIA IA IIIA IAD IIIA IA IAI AIA IAAI AIH 
IJOB = 1 
CALL FFT3D( PULSE, NCOLS , NROWS , NCOLS , NROWS ,1, IJOB, IWK, RWK, CWK) 


HAHAH IKKE KEKE EKEKKEKKEKKKKKKKKKKKKKK 


SHIFT RESULTS - IE: SWITCH QUADRANTS 1 AND 3, 2 AND 4 IN THE 
COL,ROW PLANE 


HHKKKKEKKKKKKKKKKKKKKKRKKKKKKKKKKRKKKKKKKKKKKKKRKRKKRKKKKKKKKKKKKKKRKKRKK 


CALL SHIFT (PULSE ,NROWS , NCOLS) 


VMK KHKKRKKKKEKKEKEKREKKEAKEKKERKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKkik 


MAIN PROGRAM STARTS HERE 
CALCULATIONS START AT ROW = STEP+1. 
NOTE THAT ONLY ONE QUARTER OF THE WORK MATRIX HAS TO 
BE CALCULATED. THE REMAINING THREE QUARTERS IS FORMED BY 
FOLDING AND FLIPPING. THIS IS POSSIBLE DUE TO SYMMETRY. 
HeRKRKKRKKRKRKRKRKRKRKKRKRKRKRKRKRKRKRKRKRKRKRHKRRKRAKKaaRRKRKKRRRRKRRHMKRKRRKRKRKRKHKKKKKKKRKKRK 
CALL GOTOXY(1,1) 
WRITE(*,*) *’CALCULATING ROW ’ 
IF (FTYPE .EQ. 1) THEN 
DO 200 I = STEP+1,NROWS 
CALL GOTOXY(1,18) 
WRITE(*,101) I : 
DO 210 ROW = 1,8 
DO 210 COL = 1,H 
Bl = ZPRIME(I) * RHO2(COL, ROW) 
CALL MMBSJN(B1,N,TEMP1, IER) 
WORK(COL,ROW) = TEMP1 * PULSE(COL, ROW) 
CONTINUE 


WAMAKKAKAKKKKAKKKKKKIKKAKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKAKK KKK Ik 


FLIP QUADRANT 2 TO QUADRANT 3 
FOE IOR ROR IOI TRI I TR TORT  II TOR TOR ITO TT TOTTI TOT ITT I TI TOIT TI IIIT II HH Ik 

J=0 
DO 220 ROW = H+1,NROWS 

veer) + 2 

DO 220 COL = 1,H 

WORK (COL, ROW) = WORK(COL, ROW-J) 

CONTINUE 


We Fe He He Yee He He He He He He Ve Fe Fe He He Ke te He te Fe Fe Fe He Fe He te fe te He Ve de He Fe ete Fe ic te te te te te te te ete te fete te he ete te te dete te tet tee k 
FLIP LEFT SIDE TO RIGHT SIDE 
RakKkKkKKKKKKkKKkhkkkhkhhkhhhkhkkkhkkhhkknnckknhhkhhhnhkrhehhehhhhhhhk tik 
J=0 
DO 215 COL = H+1,NCOLS 
J=uJ+2 
DO 215 ROW = 1,NROWS 
WORK(COL,ROW) = WORK(COL-J,ROW) 
CONTINUE 
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HK EKER K KKK KKK KKEAKEKCKEKKKKKKKEE 


TAKE INVERSE FFT OF EACH COLUMN. THEN TAKE INVERSE FFT OF 
ROW 33 ONLY, SINCE THIS IS THE ROW CONTAINING CENTRAL 
INFORMATION. ROW 33 BECOMES THE NEXT ROW OF THE FINAL 
OUTPUT MATRIX. 
ROTC GRICE IGISIGIGI GIS IGIGICIISIOIICIGIGIOIOIGIDI OI IGISIGIIOIEISICIDIGI GIR IGICR oIn cig iototerr 
DO 415 COL = 1,NCOLS 
DO 405 ROW = 1,NROWS 
CWK(ROW) = CONJG(WORK(COL , ROW) ) 
405 CONTINUE 
CALL FFT2C (CWK,M, IWK1) 
DO 410 ROW = 1,NROWS 
WORK (COL ,ROW) = CWK(ROW) 


ua 1 4A 4 ea eG a 


410 CONTINUE 
415 CONTINUE 
C 


DO 440 COL = 1,NCOLS 
CWK(COL) = WORK(COL,H) 
440 CONTINUE 
CALL FFT2C (CWK,M, IWK1) 


R123 NROWS * NCOLS 
C123 CMPLX(R123,0.0) 
DO 450 COL = 1,NCOLS 
RESULT(COL,I) = ABS(REAL( (CONJG(CWK(COL))/C123))) 
450 CONTINUE 


200 CONTINUE 


ELSEIF(FTYPE .EQ. 2) THEN 
TEMPL = A* Cy 2 
TEMP2 A2 * C2 / 4 
CALL GOTOXY(1,1) 
WRITE(*,*) ’CALCULATING ROW ’ 
DO 300 I = STEP+1,NROWS 
CALL GOTOXY(1,18) 
WRITE(*,101) I 
DO 310 ROW = 1,H 
DO 310 COL = 1,8 
TEMP3 = ZPRIME(I) * SQRT(ABS(RHO2(COL, ROW) **2-TEMP2) ) 
IF (ABS(RHO2(COL,ROW)) .GT. TEMP1) THEN 
CALL MMBSJN(TEMP3 ,N, TEMP4, IER) 
ELSE 
CALL MMBSIN(TEMP3,N, TEMP4, IER) 
ENDIF 
WORK (COL,ROW) = TEMP4 * EXP(-A*C2*TIME(I)) * 
& PULSE(COL , ROW) 
310 CONTINUE 
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ce, 


320 


qgqnaaa 


B15 


qg-O°0 CQ af 4 


605 


610 
615 


640 


650 


300 


FOF HWW WHIT TW HDT WHI AHHH HAHAHAHAHAHA AAAI AAA AAAI AAA AA WICK 


FLIP QUADRANT 2 TO QUADRANT 3 
ZAIDI IIR IRI III IOI III III III III III OITA TI TAA AAT AT AA HII 

J=0 
DO 320 ROW = Ht+1,NROWS 

J=J+t+2 

DO 320 COL = 1,H 

WORK (COL, ROW) = WORK(COL,ROW-J) 

CONTINUE 


Fe te He Fe te We He te We He He ee Fe He Fe He He He Fe He He Hee We TW WH KAHAN KAKKAKAHAHHK KKK AK 


FLIP LEFT SIDE TO RIGHT SIDE 
PAAR AAI TIT TI IITA IATA T IIIT TI II IIIT IIT TIT TAIT IO TIT IT T ICK 

J=0 
DO 315 COL = H+1,NCOLS 

J=Jt+2 

DO 315 ROW = 1,NROWS 

WORK (COL,ROW) = WORK(COL-J,ROW) 

CONTINUE 


te He Ye Fe He He He He He He Fe Fe Fe He He Ye He He He We Fee te Yc He He Fe He Fe He Fe Fe He Fa We He He Fe We Fe Fe Fe He Yc He Ke Fe He We He He He He ye Ye Fe He Fe He he He He He Ke Xe 


TAKE INVERSE FFT OF EACH COLUMN. THEN TAKE INVERSE FFT OF 
ROW 33 ONLY, SINCE THIS IS THE ROW CONTAINING CENTRAL 
INFORMATION. ROW 33 BECOMES THE NEXT ROW OF THE FINAL 
OUTPUT MATRIX. 
Fe Fe Fe Fe Fe Fe Fe Pe Pe BeBe RAKHI KAA AKKAAKKKAKKAHAKAKAAKKA KAA 
DO 615 COL = 1,NCOLS 
DO 605 ROW = 1,NROWS 
CWK (ROW) = CONJG(WORK(COL, ROW) ) 
CONTINUE 
CALL FFT2C (CWK,M, IWK1) 
DO 610 ROW = 1,NROWS 
WORK(COL,ROW) = CWK(ROW) 
CONTINUE 
CONTINUE 


DO 640 COL = 1,NCOLS 
CWK(COL) = WORK(COL,H) 

CONTINUE 

CALL FFT2C (CWK,M, IWK1) 


R123 = NROWS * NCOLS 
C123 = CMPLX(R123,0.0) 
DO 650 COL = 1,NCOLS 
RESULT(COL,I) = ABS(REAL( (CONJG(CWK(COL))/C123))) 
CONTINUE 


CONTINUE 


ENDIF 
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500 


© OO Qa 


510 


101 
102 


qaqaagaAAQ 


qgaQaaAa 


WH HH HHH KK KKKEKEKKEKEKEKKEKEKKEKEKEKEKEREKEEKEKKEEKEKKEKEEKEKKEKEKEKKEKKEKE 


SET ROWS 1 TO STEP TO ZERO TO SIMULATE STEP FUNCTION 
FOI IIR TOTO IOI TOR TOTTORI ROI IIIT TIT OT ITI TIT TT ITI IIR IR II 
DO 500 ROW = 1,STEP 
DO 500 COL = 1,NCOLS 
RESULT(COL,ROW) = 0.0 
CONTINUE 


We We He He He WH HHH HHH HHH WH HHH HHH HH HHH HHH KERR KKEKEEK 


OUTPUT FINAL ARRAY TO DISK 
FOTO AORTA TOR OTIC TOR ACTOR HC TOT TIO TIT TOIT TTT TT IITA TTI TTT ITI TT ITT ITI II III 
DO 510 I = 1,NROWS 
WRITE(3,102) (RESULT(J,1I),J = 1,NCOLS) 
CONTINUE 


FORMAT (12) 
FORMAT (64F16.7) 


te te ve He te Ye Ve ve He He We ve He ve He He ve ve He ve He ve Ye He ve ¥e Ye We ve Ae oe He He He Ye Ye Ye Ye ve He He He ve Ye Ye He He te Ye Ye Ve Fe Ye He He He He ve ¥e Ye ve ve He ee 


CLOSE REMAINING OPEN FILES 


te te ve Ye ve Fe te He He ve te He ve He He He 40 ee 0 He He He He He ie ie Ye 0 0 0 0 0 0 ie Hee Yee ie et ok ek ie ek ete tere ek te ke Wee He ok 


CLOSE(2) 
CLOSE(3) 


He Ve He He We He He Wee He He He HH HH He IK HH KK HHH HIRI KIKI HHH K KEK RAK 


REMIND USER OF DATA FILES CREATED 

Ye ve te Fe Yo Ye ve Ve He He te He ve He We eve He Xe He ¥e Ye He He Fe He He te te He 00 Fe 0 ee ee HH KK He te RAK KK KITE 
CALL CLRSCR 

WRITE(*,*) ’FINISHED.’ 

WRITE(*, *) 

WRITE(*,*) ’ THE FOLLOWING ASCII TEXT FILES ARE IN THE DEFAULT’ 
WRITE(*,*) ’ DIRECTORY: ’ 

WRITE(*, *) 

WRITE(*,*) ' FILENAME DESCRIPTION’ 

WRITE(*,*) § 0 --------+------------------------------------------ 
WRITE(*,*) ’ ',FNAME1,’ INPUT IMPULSE FUNCTION’ 
WRITE(*,*) ’ ’,FNAME2,’ SUMMARY INFORMATION’ 

WRITE(*,*) ’ ’,FNAME3,’ OUTPUT ARRAY’ 

WRITE(*, *) 

WRITE(*, *) 


STOP 
END 
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He te HK HH HK KKK KHER KEKEKKKEKKKEKKEKKEKEKKEK 
* * 
* SUBROUTINE GOTOXY 
* MOVES CURSOR TO LINE X, COLUMN Y * 
CALL: CALL GOTOXY(X,Y) * 


* * 


He He KHAKI KKH KE KKK EEE KEKEEKEKEKEKE KEKE KEKKEKKKKEKKKKEK 


qaQaaaaaaa 


SUBROUTINE GOTOXY(X,Y) 

CHARACTER*1 C1,C2,C5,C8,LC(5) 

CHARACTER*5 CBUFF 

INTEGER*2 IC(4),L,X,Y 

HOUIVALENGE (C1,1C(1)),(C2,1C(2)),(C5,1C(3)),(C8,ICC4)), 
SCCBhUPE,LCO(C1)) 

DATA IC/16#1B, 16#5B, 16#3B, 16#66/ 


L = 10000+100*X+Y 


QD 


*** WRITE ESCAPE CODES TO CHARACTER BUFFER 
WRITE(CBUFF,2) L 
FORMAT(I5) 


aA & 


QO 


*** WRITE ESCAPE CODES TO THE DISPLAY 
WRITE(*,1) €C1,C2,LC(2),LC(3),C5,LC(4),LC(5),C8 
-FORMAT( 1X, 8A1, \) e 


me 


Q 


*** END OF SUBROUTINE 
RETURN 
END 


Ke HM He He HHH KH KH KKK KKK KK KKK HK KKK KEKE KKKKKKEKKKKKKKKKKKKKKEKEK 
* ve 
* SUBROUTINE CLRSCR * 
* SUBROUTINE TO CLEAR THE DISPLAY * 
* CALL: CALL CLRSCR * 
ve * 


REMAKE KKKKEKKKKEAKKKEEKEHKKEEKKEKAE KEKE KKK KKKEKK KKK KKKKKK 


qgaaQaaaaaa a 


SUBROUTINE CLRSCR 

CHARACTER*1 C1,C2,C3,C4 

INTEGER*2 IC(4) 

EQUIVALENCE (C1,1C(1)),(C2,1C(2)), (C3,1C(3)), (C4, IC(4)) 
DATA IC/16#1B, 16#5B, 16#32, 16#4A/ 


GC *** WRITE ESCAPE CODE TO DISPLAY 
WRITEC*, 1) C1,CZ,C3,C4 


i FORMAT (1X, 4A1) 

C 

C *** END OF SUBROUTINE 
RETURN 
END 
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qeanNgangaNa|gndaaa aA 4a 


QO 


10 


FFF HH KKK III HIKER KKH RKKKKKKKKKKKK 


* SUBROUTINE TO SHIFT A SQUARE MATRIX - EXCHANGES QUADRANTS ONE * 


* AND THREE, AND TWO AND FOUR. 
* USE: CALL SHIFT(X,NROWS ,NCOLS) 


: WHERE X IS A SQUARE ARRAY - REAL OR COMPLEX 
= NROWS ,NCOLS IS THE ARRAY SIZE 


* 


* 


* 


* 


* 


* 


Fe He He Fe Ve He Fe Ve Ve Ve Ve Ve Ve Ve HI HHI KKH KKK KKK KKK KKK KEKE HEHEHE KKKKK 


SUBROUTINE SHIFT (X,NROWS,NCOLS) 


INTEGER NROWS ,NCOLS , ROW,COL,R2Z,C2 
COMPLEX X(NCOLS,NROWS) ,T1,T2 


R2 
C2 


NROWS / 2 
NCOLS/2 


DO 10 ROW = 1,R2 
DO 10 COL = 1,C2 

Ti = X(COL, ROW) 
X(COL,ROW) = X(COL+C2, ROW+R2) 
X(COL+C2Z ,ROW+R2Z) = T1 
T2 = X(COLtC2, ROW) 
X(COL+C2,ROW) = X(COL,ROW+R2) 
X(COL,ROW+R2) = T2 

CONTINUE 

RETURN 

END 
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